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Abstract. We have developed a method based on the embedded Kadanoff-Baym equations 
to study the time evolution of open and inhomogeneous systems. The equation of motion for 
the Green's function on the Keldysh contour is solved using different conserving many-body 
approximations for the self-energy. Our formulation incorporates basic conservation laws, such 
as particle conservation, and includes both initial correlations and initial embedding effects, 
without restrictions on the time-dependence of the external driving field. We present results 
for the time-dependent density, current and dipole moment for a correlated tight binding chain 
connected to one-dimensional non-interacting leads exposed to DC and AC biases of various 
forms. We find that the self-consistent 2B and GW approximations are in extremely good 
agreement with each other at all times, for the long-range interactions that we consider. In the 
DC case we show that the oscillations in the transients can be understood from interchain and 
lead-chain transitions in the system and find that the dominant frequency corresponds to the 
HOMO-LUMO transition of the central wire. For AC biases with odd inversion symmetry odd 
harmonics to high harmonic order in the driving frequency are observed in the dipole moment, 
whereas for asymmetric applied bias also even harmonics have considerable intensity. In both 
cases we find that the HOMO-LUMO transition strongly mixes with the harmonics leading to 
harmonic peaks with enhanced intensity at the HOMO-LUMO transition energy. 



1. Introduction 

The progress in miniaturization of electronic devices has led to a completely new field of research, 
known as molecular electronics, in which it has become possible to prepare samples and to 
measure steady state conductivity properties of systems with single molecules squeezed between 
conducting electrodes [U [2]. However, in future applications of nanodevices one is not only 
interested in their steady properties, but one also wants to manipulate these systems in time. For 
the control and operation of such devices the understanding of the time-dependent properties, 
such as switch-on times, currents, density fluctuations etc., becomes crucially important. In 
contrast to the steady-state approach, understanding these phenomena demands truly real-time 
treatment of the problem which is the topic of this paper. 

The study of periodic pulses is of special importance to manipulate and control the electron 
flow in nanoscale devices. In this paper we present results obtained from the solution of the 
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Figure 1. A) The transport setup under consideration: A correlated central region is coupled 
to semi-infinite left and right noninteracting tight-binding leads. B) The Keldysh contour 7. 
Times on the upper /lower real-time branches are denoted with t^. The imaginary (thermal) 
track runs from Iq to to ~ with f3 the inverse temperature. 

Kadanoff-Baym (KB) equations for periodic voltage pulses. Within this framework also periodic 
gate voltages or barriers can be included without any additional computational effort. All these 
kind of signals can be experimentally realized and have been used to pump the electron current in 
several nanoscale structures ranging from molecules and quantum dots [3lll] to nanotubes [5j 6j. 
The AC transport properties of nanoscale systems have been discussed recently [71 |Hl El |10] and 
moreover, the advantages of real-time approaches in the context of quantum transport through 
molecular junctions driven out of equilibrium by periodic fields have been discussed in detail in 



There is another reason for studying AC fields in quantum transport. Usually the assump- 
tion is made that the external field in the leads is instantaneously screened. However, in time- 
dependent transport transient times can be of the same order as the plasma oscillation period 
(see e.g. Ref. [8J). In such a case it would be reasonable to assume that at Hartree mean- field 
level the bias in the leads can be described by an AC field that alternates with the plasma 
frequency. 

In this article we use a recently proposed theoretical framework \12\ [T3] to study time- 
dependent electron transport through a quantum wire connected to two one-dimensional metallic 
leads in the presence AC and DC fields. The theoretical approach is based on the real-time 
propagation of the KB equations [HI [151 HSl [13 [El [IHl 120] for open and interacting systems. 
The electron-electron interactions are handled perturbatively through the many-body self-energy 
for which we have implemented the HF (Hartree Fock), 2B (Second Born) and GW conserving 
approximations while the leads are taken into account through embedding self-energies. We 
can calculate a wealth of observables such as currents, densities and dipole moments, and study 
electronic and initial correlations induced by the time nonlocal self-energy. 

2. Theory 

2.1. The Hamiltonian 

We consider a quantum correlated open system (a central region) coupled to noninteracting 
electron reservoirs (leads) and described by the Hamiltonian (see Fig[TjA) 
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Here He, Ha and Ht are the central region, the lead a = L,R and the tunneling Hamiltonians 
respectively, and N is the particle number operator coupled to the chemical potential fi. The 



explicit expressions for these Hamiltonians read 



ij,a ijkl 

Kit) = U^{t)Na + Y,h'^clj,^o., (3) 

Hj, = ^ Vija[dl^Cj„a + c]^Jia], (4) 
2 J, era 

where a, a' are the spin indices, d^d and c^,c are the creation and annihilation operators in the 
device and lead region and i,j label a complete set of one-particle states in the corresponding 
subspaces of the system. The one-body part hij{t) of the central region Hamiltonian is in 
general time-dependent to account for driving fields such as gate voltages or pumping fields. 
The terms vijki in the two-body part are the two-electron integrals of the Coulomb interaction. 
Furthermore, the one-body part hf- of the lead Hamiltonian describes the metallic leads and 

the tunneling Hamiltonian Ht contains the lead couplings to the central region. The system is 
driven out of equilibrium by a homogeneous time-dependent bias voltage Ua{t) coupled to the 
number operator = Yli a ^La^*""" for lead a. This field can be considered as the sum of the 
externally applied field and the resulting screening field of the metallic leads which could, for 
instance, represent an AC field describing a plasmon oscillation. 



2.2. Kadanoff-Baym equations for the Keldysh Green's function 

We study a system that is in equilibrium at an inverse temperature /3 at times t < to and 
described by a Hamiltonian Hq. For t > tQ the system is driven out of equilibrium by an 
external bias voltage and evolves in time according to the Hamiltonian H{t). The Keldysh 
Green's function for this system is defined as the expectation value of the contour-ordered 
product of the creation and annihilation operators [151 IISl ED [22] 

Qrsiz, z') = -i ^ ^, (5) 

where a and a) are either lead or central region operators and the indices r and s are collective 
indices for position and spin. Furthermore, z is a contour time variable and T orders the 
operators along the Keldysh contour (see Fig{T]B) by arranging the operators with later contour 
times to the left. The trace is taken over the many-body states of the system. All the time- 
dependent one-particle properties can be calculated from the Green's function, for example the 
time-dependent density matrix is given by 

nrs{t) = -iQrs{t-,t+), (6) 

where the time arguments t± are located on the lower/upper branch of the Keldysh contour. 
Starting from Eq. (0) the equation of motion for the Green's function and its adjoint can easily 
be derived and are given by 

id,G{z,z') = d{z,z')i + n{z)g{z,z')+ [ dz-s^^{z,z)g{z,z'), (7) 

-id,>g{z,z') = 6{z,z')l + g{z,z')U{z')+ [ dzg{z,z)'S^^{z,z'), (8) 



where S'^^ is the many-body self-energy, H(2) is the one-body part of the fuh Hamiltonian and 
the integration is performed over the Keldysh-contour. These equations of motion need to be 
solved with the Kubo-Martin-Schwinger (KMS) boundary conditions: G{to, z') = —G{to — il3, z') 
and Q{z^ to) = —Q{z, to—i/S). All the quantities are considered as block matrices with a structure 
given by 
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(9) 

where the different block matrices describe the projections onto different subregions. We assume 
that there is no direct coupling between the leads. The many-body self-energy has nonzero 
elements only in the central region. This follows immediately from the diagrammatic expansion 
of the self-energy. More specifically, if we expand the self-energy in powers of the two-body 
interactions we see that the interaction matrix elements Vjjfc; only connect sites in the central 
region and therefore all Green function lines in the corresponding Feynman diagrams are of 
the type Qqq- Therefore 5]^^[^cc] is a functional of Qqq only. Note that we do not expand 
in the lead-device couplings which are all exactly incorporated in the one-body matrix H. To 
study the dynamical processes and to extract the observables of interest such as the densities 
and currents in the whole open and connected system we need the solution for the full Green 
function in different subregions of the system. In particular we need the Green's functions Qqq 
projected onto central region C and Qaa projected onto the lead a region. For these we need to 
extract from the block matrix structure of the Green's function an equation for Qqq and Qaa. 
The projection of the equation of motion ([7]) onto regions CC and aC gives 

[ld,\-YLcG{z)\GGQ{z,z')=b{z,2!)\ +^HcaC?ac(^,^')+ / dzYl^^^iz ,z)G gq{z, , 

a. h 

(10) 



for the central region and 



jia^l - Yiaa{z)^Qac{,Z,z') = HacG Cc{z , z') , 



(11) 



for the projection on the interface aC region. We see that the right hand side of Eq. (|lip is of a 
simple form as a result of the absence of two-body interactions in the leads. The uncontacted 
and biased lead Green's function g^^ satisfies the equation of motion 



^id^l - liaaiz)^gaaiz, z) = 6{z, z')l. 

Using this equation we can construct the solution for GaC and it reads 

Qac{z,z')= / dzg^a{z,z)UacQcc{z,z'). 

J-y 

Taking into account Eq. (|13|) the second term on the righthand side of Eq. (jlO|) becomes 



dz,z 



/ dzT,cm{z,z)gcc{z,z'), 



(12) 



(13) 
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where we have introduced the embedding self-energy 



^cm{z,z') = '^'Ecm,aiz,z') = '^Uca 9 aaiz , z')UaC , 
a a 



(15) 



accounting for the tunneling of electrons between the central region and the leads. As we see 
from the definition the embedding self-energy Sem,a depends only on the coupling Hamiltonians 
and on the isolated lead Green's function g^a- Furthermore, the Green's function is determined 
once the isolated lead Hamiltonian Ha of Eq. ([3]) is specified. Inserting ()14p back to (jlOp then 
gives the equation of motion 

{id,l-Ucc{z)}Gcc{z,z') = 5{z,z')l + J^dz + {z,z)Gcc{z,z'). (16) 

An adjoint equation can be derived similarly [13] . Equation (jl6p is an exact equation for the 
Green's function Gcc for the class of Hamiltonians of Eq. ([1]), provided that an exact expression 
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as a functional of Qcc is inserted. 



To apply Eq. ()16p in practice we need to transform it to real-time equations that we solve by 
time-propagation. This can be done in Eq. ()16p by considering time-arguments of the Green's 
function and self-energy on different branches of the contour. We therefore have to define these 
components first. Let us thus consider a function on the Keldysh contour of the general form 



F{z, z') = F\z)6{z, z') + e{z, z')F>{z, z') + e{z\ z)F<{z, z'), 



(17) 



where 6{z,z') is a contour Heaviside function [19] . i.e., 0{z,z') = 1 for z later than z' on the 
contour and zero otherwise, and 6{z, z') = dz0{z, z') is the contour delta function. By restricting 
the variables z and z' on different branches of the contour we can define the various components 
of F as 



F^it,t') 

F^t,T) 

F^ir,t) 
F'^T-r') 



F{t±,to - it), 
F{to - IT, t±), 
-iF{to - IT, to - it'), 



(18) 
(19) 
(20) 
(21) 



and 

F^M(t, t') = F\t)5{t - t') ±e{±tT t') [F> {t, t') - F< (t, t')] . (22) 

For the Green's function there is no singular contribution, i.e., = 0, but the self-energy has 
a singular contribution of Hartree-Fock form, i.e., 'S^ = S^^[^] [19j. With these definitions 
we can now convert Eq. ()16p into equations for the separate components. This is conveniently 
done using the conversion table in Ref . |23j . This leads to a set of coupled real-time differential 
equations, known as the Kadanoff-Baym equations 



idtg^{t,t') 
-a.C?^^(r-T') 



g^{t, t')Hcc(t') + Ig"^ • + c?^ • + {t, t'), 



C?r(r,t)Hcc(t)+ + (T,t), 

15(r - r') + Uccg^'^ir - t) + i [E^'^ * g''] (r - r'). 



(23) 
(24) 
(25) 
(26) 
(27) 



Here the self-energy is the sum of the many-body and embedding self-energies. The superscript 
notation for the various components of self-energy and the Green's function, "M",<,>,] and [ 
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Figure 2. Diagrammatic representation of the self-energies for the conserving many-body 
approximations used. 

are used to identify the time-arguments on different parts of the Keldysh contour [1'6\ \T2\ [23] . 
The notations • and * denote the real-time and imaginary-time convolutions: 

/•OO 1-0 

[a-b]{t,t')= / a{t,t}b(t,t')dt , [a-kb]{t,t') = -i j a{t,T)b{T,t')dT. (28) 
^0 Jo 

In practice, the Matsubara equation (|27p is solved first since it is disconnected from the real- 
time KB equations. The real-time Green's functions are then initialized with the Matsubara 
Green's function as: C?>(0,0) = iC?^(0+), C?<(0,0) = iG^^Q-), G\o,t) = iQ^^-t) and 
Q^{t,Q) = iQ^^ [t) and then the equations are solved using the time-propagation 

method described in Ref. 



2.3. Conserving approximations 

The importance of particle number conservation in quantum transport transport has been 
carefully addressed in Refs. I26j . The particle number conservation law in open systems 
is given by 

^ = /L(t) + lR(t), (29) 

and expresses the fact that the change in the number of particles in the central region sums up 
to the total current that flows into the leads. Particle number conservation is guaranteed within 
the KB approach by using expressions for the self-energy that can be obtained as functional 
derivative with respect to the Green's function of a functional <I>: 

5^CC,rs [yCC\ [Z,Z) = — — — r . (30) 

Baym has proven ^7] that such ^-derivable or conserving approximations |28t [271 ES] 
automatically lead to satisfaction of the conservation laws provided that the equations of 
motion for the Green's function are solved to full self-consistency. Commonly used conserving 
approximations for I]q,?[^cc] are the Hartree-Fock, second Born and GW approximations 
[151 [Ml [SQl [31] which are displayed diagrammatically in Fig[2l 



2.4. Time- dependent current 

An equation for the time-dependent current flowing into the lead a can be derived from the 
time-derivative of the number of particles in lead a using the equation of motion for the lead 



Green's function Q^a [IS]- This yields 



= 2ReTrc [g-^^{t,t)H^c\ 



(31) 



If we insert the adjoint of Eq. (|13p and extract the different components from the resulting 
equation using the conversion table in Ref. |23j . the equation for the current becomes [13] 



I^{t) = 2ReTrc 



J to J to 



[ (iT4c(*'^)^L,a(T,t) 
JO 



(32) 



The first two terms in this equation contain integrations over earlier times from to to t and take 
into account the nontrivial memory effects arising from the time-nonlocality of the embedding 
self-energy and the Green's functions. Furthermore, the last integral over the imaginary 
track incorporates the initial many-body and embedding effects. Equation (|32|) provides a 
generalization of the widely used Meir-Wingreen formula |32[ [55] to the transient time-domain 
ll3j. 



2.5. Nonequilibrium spectral function 

The nonequilibrium spectral function for real times is defined as 

Ait,t')=iTTc[G^c-Gcc]M= E i{diAt),dlAt')]). (33) 

where the brackets in the last term denote the anti-commutator. It is common to Fourier 
transform this function with respect to the relative time-coordinate t = t — t' for a given value 
of the average time T = (t + t')/2. The spectral function in the new coordinates is then 

^(r,a;) = -Im J ^ e'-Trc[0gc " ^ccK^ + ^> ^ " ^) = Re J ^e'^^ A{T +^-,T -^-) (34) 

If the time-dependent external field becomes constant after some switching time, then also the 
spectral function becomes independent of T after some transient period. In that case we will 
denote the spectral function by A{ijj). The spectral function for a system in equilibrium has 
peaks at the addition and removal energies. In the nonequilibrium case this is less clear as one 
can, in general, not use a Lehmann expansion to prove this fact. However, for a system to which 
a bias is suddenly applied one can still show that the spectral function will contain peaks at the 
addition and removal energies from the various excited states of the biased system. 



2.6. The density in the leads 

The density in the leads can be obtained from the Green's function Qaa- This Green's function 
satisfies an equation of motion that can be obtained by projecting the equation dTj) into the aa 
region: 

^id:,l - li.aa{z)^Qaa{z,z') = 5{z,z')l H„c^Ca (^, ^0 • (35) 

Using the adjoint of (fT5|) in ([5^ we have 

^idz -'H-aa{z)^Qaa{z,z') = 5{z, z')l + j dzYli^^^iz , z)g ^^z , z') , (36) 



where we defined the inbedding self-energy as 



'^m,ai^, z') = UacGcciz, z')'H.Ca- (37) 

Equation ()36p can now be solved using Eq. (|12p . By taking the time arguments at t± we get 

Gaait-,t+) = g^^{t-,t+) + dzdzg^^{t^,z)'Sir,^o.{z,z)g^^{z,t+). (38) 

We see from Eq.© that from this equation we can extract the time-dependent spin occupation 
of orbital i in lead a by taking r = s = iaa. The first term on the r.h.s. of Eq. ()38p gives 
the density of uncontacted biased lead a and the second term gives a correction to the den- 
sity induced by the presence of the correlated scattering region. In practice we first solve the 
Kadanoff-Baym equations for Qcc and then we construct the inbedding self-energy which 
is then used to calculate the time-dependent lead density according to (f38|) . 



2.7. Embedding self-energy 

In this section we derive the expression for the embedding self-energy for the case of one- 
dimensional noninter acting tight-binding leads. Prom equation ()15p we see that the embedding 
self-energy is given by 

Afc, 

{z, z') = Vk ,ia9aa,ij{Zi Z )Vja,l, (39) 

where k,l label the sites in the central region and i,j the sites in the lead a. Furthermore, J\fa 
is the number of sites in the lead a (at the end of the derivation we take Afa — > oo). We now 
express the biased lead Green's functions using the eigenbasis of the biased lead Hamiltonian: 

pq 

where 

~9tc.,pq{z,z') = ^5p,/(ep„)e-/.'(^--'^+^'^(^"))'^^ (41) 
9^a,pqiz,^') = ^'5p,(/(6pa)-l)e-*/^'(^--^+^'^(^"))'^^ (42) 

are the lead Green's functions in the eigenbasis of the leads with eigenvalues Cpa- The columns 
of matrix U" represent these eigenstates in site basis. The function /(e) = l/(e^('^~^'^ -|- 1) is the 
Fermi distribution function. Inserting Eqs. (j4ip and (j42p into the definition of the embedding 
self-energy ()39p and defining the function 

TmA^) = 2^ ^kMK^i^ - epa)^^;/ (43) 

ijp 

we get an expression for the one-dimensional lead self-energy 



The remaining task is to evaluate the function Tki^ai^)- For a tight-binding (TB) chain we have 



a" + 2raCos((/)p), (46) 



where a" and Tq, are the tight binding on-site and hopping parameters for the lead a and 
~ aTTT' P ~ ^■■■■^a = I.-.A/'q. We consider the case that the coupling Hamiltonian has 
nonzero elements only at the terminal points adjacent to the endpoints of the central region 
nanowire, that is when i = j = 1. The function Tki^ai^) then becomes 

TmA^) = '^^Vk,iaVio.,i 5^ f/iV(e - ep^)U;l (48) 

p=i 

Inserting Eqs. ([I7|) into ([l8]l and taking the Ma — > oo limit we obtain 



rMA^)= j£| y^-^^J x0(2|r,|-|e-a"|). (49) 

From the definition of rfc/ „(e) we see that it describes the weighted density of states of the 
lead a as a function of the energy. Equation (j49p then yields the expression for the ID lead 
self-energy 



7r|Ja| Ja'--2\Tc\ V V ■^^a J 

and the ^ component is obtained by simply replacing the /(e) with /(e) — 1. Furthermore, 
the [,] and "M" components are similarly obtained by considering the time-arguments on 
different parts of the Keldysh contour. 

3. Numerical results 

We now apply the KB approach to study the transport dynamics of a quantum wire consisting of 
4 sites connected to left and right one-dimensional TB leads, see Fig. [T]A. The nearest neighbour 
hoppings in the leads are set to Tl = Tr = — 2.0 and the on-site energy a" is set equal to the 
chemical potential, = /u, yielding a half-filling for the lead energy bands. The values for 
the coupling Hamiltonian are set to Vi^il = V^^ir = —0.5 and the central region sites have the 
on-site energies ha = and hopping parameters hij = —1.0 between nearest neighbour sites i 
and j. The electron-electron interaction in the central region has the form vij^i = Vij SuSjk with 



where vu = 1.5. The chemical potential is fixed between the highest occupied molecular orbital 
(HOMO) and lowest unoccupied molecular orbital (LUMO) levels of the isolated central region 
and has the value fi = 2.26. Furthermore, the inverse temperature is set to /3 = 90 corresponding 
to the zero temperature limit. We use different time-dependent bias voltages. We consider a 
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Figure 3. Left panel: Time-dependent right current I(t) and dipole moment d{t) for different 
self-energy approximations when a DC bias Uq = 1.4 is applied to the leads. The inset shows 
the absolute value Fourier transform of the dipole moment. Right panel: Steady state spectral 
function A{u}). 
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Figure 4. Panel a): Dipole moment d{t) corresponding to bias Ui{t). The inset shows the 
absolute value of the Fourier transform of the dipole moment. Panel b): Time-dependent current 
flowing into the right lead (top panel) and the time-dependent bias Ui{t) (bottom panel). 



sudden switch on DC case of the form UL{t) = —UR^t) = Uo6{t — to) and periodic AC pulses 
Uiit) = -Unit) = ?7i,2(t) of the forms 

Ui{t) = (block bias) (52) 




<t<Tr 



0<t<^ 



U2{t) = { ^ ^- 2 (on-offbias) (53) 



which for t > Tq are periodically extended using U{t + Tq) = U{t). The time Tq is the period 
of the complete bias cycle corresponding to oscillation frequency loq = 2-k/Tq. In all cases we 
choose the amplitude of the bias voltage to be C/q = 1.4 and for AC biases we set the period to 
be Tq = 40 a.u. Finally, to study how the charge distributes along the chain after a bias voltage 




Figure 5. Panel a): Dipole moment d{t) corresponding to bias U2{t). The inset shows the 
absolute value of the Fourier transform of the dipole moment. Panel b): Time-dependent current 
flowing into the right lead (top panel) and the time-dependent bias U2{t) (bottom panel). 



is switched on we calculate the time-dependent dipole moment 

4 

d{t) =^Xini{t) (54) 
1=1 

where the Xj are the coordinates of the sites of the chain (with a lattice spacing of one), with 
origin halfway between sites 2 and 3, and nj(t) is the site occupation number. 
We start by considering the case of a suddenly switched DC bias. In the left panel of Fig. [3] we 
compare the time-dependent dipole moments (three lowermost curves) and transient currents 
(three uppermost curves) for the HF, 2B and GW approximations. The right panel shows the 
steady-state spectral functions. We see that the Hartree-Fock peaks are relatively sharp with 
a width approximately given by the imaginary part of the embedding self-energy at the Fermi 
level of r = 21/ ~ 0.25 where V = Vi^n = Va^ir. On the other hand the 2B and GW spec- 
tral functions are very broadened due to strongly enhanced quasi-particle scattering at finite 
bias [26t JSj . We see that the Hartree-Fock approximation produces a slightly larger current 
compared to the 2B and GW approximations. This is in keeping with the fact that the HF 
spectral function integrated over the bias window leads to a slightly larger value as compared 
to those from the 2B and GW approximations. We further see that the correlated 2B and GW 
approximations lead to currents and dipole moments in close agreement with each other. This 
indicates that the dynamical screening of the Coulomb interaction as described by the first bub- 
ble diagram of the self-energy has a leading contribution. 

In the inset of the left panel in Fig. [3] we display the absolute value of the Fourier trans- 
form of the dipole moment. We see a number of peaks that can be directly related to transitions 
from the lead electrochemical potentials fi + Ua to the lead levels and between levels in the 
central sites themselves. The nature of these transitions can be deduced from the peak positions 
in the spectral function. Denoting these energies with Aen, Aej/j and Aij we can identify, at HF 
level, transitions at energies A.eL3,2R ~ 0.5, AeL4 0.9, A23 ~ 1.6, A34 « 1.3, AeL2,3_R ~ 2.25 
and AeLi,4R ~ 3.3. The largest peaks are due to the HOMO-LUMO transition A23 and the 
transitions Aen, AeL3,2R and we see that these are the only peaks that are clearly visible for 
the 2B/GW case, since in the correlated case the transients and dipole moments are strongly 
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Figure 6. Time-dependent spectral function A(T, u) in the HF (left panel) and 2B (right panel) 
approximation. The system parameters are the same as for Fig. [5j 

damped. 

We now address the case of the AC block bias Ui{t). In Fig. SI we show the corresponding 
dipole moments and transient currents. As in the DC bias case we observe that the GW results 
are very close to the 2B results. The currents and dipole moments show a similar initial tran- 
sient behavior as in Fig. [3] up to t = 20 after which the bias switches sign. The inset in the 
Fig. 0] a) contains the Fourier transform of the dipole moment which displays peaks at the odd 
harmonics (2n -|- l)u;o of the basic driving frequency ujq of the AC bias. The appearance of odd 
harmonics is a consequence of the even inversion symmetry of the unbiased Hamiltonian and 
the odd inversion symmetry of the applied bias. We see up to 9 harmonic peaks which indicates 
that we are far from the linear response regime. From the inset of Fig. H] we also see that the 
harmonics at energies 1.4 (harmonic order 9) and 1.7 (harmonic order 11) are enhanced. This 
is due to strong mixing with the HOMO-LUMO transition that occurs at energy 1.6 in the DC 
case (see figED- 

High harmonics with even order will develop when we break the symmetry of the unbiased 
Hamiltonian or of the applied bias. This happens for the case of the on-off AC bias voltage 
U2(t) (see Eg. ([55]) ). The dipole moments and transients for this case are displayed in Fig. EJ 
The applied bias in this case has either even or odd inversion symmetry depending on the time 
interval considered. Therefore even harmonics at frequencies 2noJo are also visible in the Fourier 
transform of the dipole moment (Fig. [5] inset). Also in this case the higher harmonics mix 
strongly with the main electronic transition frequencies leading to an increased intensity of the 
harmonics around the HOMO-LUMO gap frequency u! ~ 1.5. Furthermore there are increased 
harmonic peaks at frequencies uj ~ 0.5,0.9 which are due to mixing of the harmonics with the 
transitions AeL3,2R and AeLi (see figl3]). 

Now we turn our attention to the spectral functions. In Fig. [6] we show the time-dependent 
spectral function A{T,uj) (per spin and shifted with the chemical potential fi) for the HF and 
2B approximations for the case of the on-off bias U2{t). The system is first propagated without 
the external bias up to 40 a.u. after which the AC bias is turned on. We see that the spectral 
peaks in Hartree-Fock approximation remain sharp during the whole time-propagation and due 
to the applied bias a closing of the HOMO-LUMO gap takes place. In contrast to the HF ap- 
proximation, the spectral peaks in 2B start to broaden and loose intensity rapidly after the bias 
has been switched on and is caused by enhanced quasiparticle scattering at finite bias. We note 
that we do not observe any noticeable shift in the position of the spectral peaks as a function of 




Figure 7. Left panel: Real-time spectral function A(t, t') in the 2B approximation for on- 
off AC bias voltage. Right panel: The real and imaginary parts of A(T + t/2,T — t/2) as a 
function of r/2 for T=70 (solid line) and T=90 (dashed line). The corresponding cross-sections 
are indicated with black lines. 



the alternating bias. This can be explained from an analysis of the spectral function A{t,t') in 
real time. We display contour plots of these functions in Figs. [7] and [SI In these figures we also 
display two center-of-time cross sections A{T + t/2,T — t/2) for T = 70 and T = 90. Also for 
this case the ground state was first propagated without a bias up to 40 a.u. after which the AC 
field was switched on. Therefore for t,t' < 40 the contour plot describes the equilibrium system. 
For such times the spectral function only depends on the relative time t — t'. For t, t' > 40 the 
bias is applied and the spectral function starts to depend on both time-coordinates separately. 
This is especially visible for the 2B case in FiglT] where for 40 < t, t' < 50 we see a sudden 
collapse of the spectral function (we have checked this also for other values of the contour lines 
than displayed in the plot). This time interval coincides with the transient time of the current as 
shown in Figl5l Note, however, that the spectral function A(T,uj) for average times 20 < T < 30 
is still dominated by the equilibrium state since the Fourier transform over the relative time r 
extends mainly over values of t, t' within the equilibrium region of the double time plane. For 
this reason the collapse of the spectral function A{T,uj) of Figl6] occurs within the time interval 
from T = 20 to T = 40. In Fig. [7] we further observe a characteristic checkerboard pattern of 
the contour line at zero. This reflects the various combinations (on-on, on-off and off-off) of the 
functions U2{t) and U2{t') at times t and t' . If we calculate the spectral function A{T,uj) by 
Fourier transforming over the relative time coordinate r we make an average over these regions. 
This averaging makes the spectral function less sensitive to the period of the applied AC bias. 
Nevertheless, some periodic changes in the cross sections A{T + t/2, T — t/2) can be observed. 
In the upper and lower right panels of FiglT] we display these cross sections at average times 
T = 70 (center of the off-period) and 90 (center of the on-period) for both the real and the 
imaginary part of A(T + r/2,T — t/2) as a function of r/2. On the time diagonal the real 
part of A is always equal to 4 since we have four states per spin in our atomic chain. Since 




Figure 8. Left panel: Real-time spectral function A{t, t') in the HF approximation for on-off 
AC bias voltage . Right panel: The real and imaginary parts of A{T + t/2,T — t/2) as a 
function of r/2 for T=70 (solid line) and T=90 (dashed line). The corresponding cross-sections 
are indicated with black lines. 



the spectral function decays fast as a function of r the main changes in the spectral function 
occur close to the time diagonal. When we Fourier transform these cross section small shifts 
in the spectral peaks are observed. In the on-state at T = 90 we found a slight closing of the 
HOMO-LUMO gap as compared to the off-state at T = 70. 

In FiglS] we display the spectral function A(t, t') for the HF approximation. It can be seen 
that the oscillations along the relative time coordinate r in the spectral functions are much less 
damped in HF (Fig. [8]) than in 2B (Fig. [7]). As a consequence we do not see the checkerboard 
pattern observed in the 2B case since it occurs much farther away from the time diagonal. The 
functions A{T + r/2,T — t/2) clearly show a slow and a fast oscillation as a function of t/2 
corresponding to the two outer and inner peaks (at positive and negative frequencies) of the 
time-dependent spectral function A{T,uj) displayed in Fig. [6l Note that in this figure we shifted 
the peaks with the chemical potential fj, = 2.26. 

We now turn our attention to the lead densities for the case of the AC block bias Ui{t). In 
right/left panel of Fig. [9] we show the time evolution of the site occupation numbers ni{t) in the 
right lead for the AC/DC biases. These are calculated for up to 30 sites deep into the lead from 
Eq. ([38|) using ni{t) = —iQaa,rr{t--,t+) where r = ia. For the case of a DC bias shown in panel 
(a) we can clearly observe a density wave front that propagates into the right lead. Furthermore 
the electron density displays a zigzag-pattern with an amplitude that decreases away from the 
first site into the lead. These are the usual Friedel oscillations caused by the presence of the 
atomic chain in the lattice. Their spatial profile is given by 5n{x) ^ sm.{2kpx) where kp is the 
Fermi wave vector and x the distance to the impurity site. In the case of the half-filled lead 
energy band that we consider kp = vr/2 corresponding to a density variation that alternates on 
every site. When we switch on the bias in the leads the Friedel oscillations get enhanced by 




Figure 9. Site occupation number nj(t) per spin in the right lead as a function of time t and 
site label i. The results are displayed for (a) DC bias, and (b) AC bias Ui{t) within the HF 
approximation. 

roughly an order of magnitude. 

In panel (b) of Figi9]we shown the lead density profile for the case of the AC bias Ui{t). As in 
the DC case, we see a wavefront moving into the lead. However, due to the periodic alternation 
of the bias these wavefronts get superimposed on the returning ones giving rise to interferences 
causing additional wiggles in the density profile. 

4. Conclusions 

We proposed a time-dependent many-body approach based on the real time propagation of the 
KB equations for open and inhomogeneous systems. 

The method was applied to study transport dynamics through an atomic chain under AC 
and DC bias voltages. We calculated the spectral functions and the time-dependent current, 
density and dipole moment within the HF, second Born and GW conserving approximations. 
We found that electron correlations beyond mean field have a large impact on time-dependent 
and steady state properties. Both in the AC and DC case the HF spectral function retains its 
sharp structures upon applying a bias, while the spectral function calculated within the 2B and 
GW approximations broadens considerably. 

The strong AC biases that we applied lead to highly nonlinear perturbations. As a 
consequence high-order harmonics of the driving frequency are observed in the time-dependent 
dipole moment. Depending on the symmetry of the Hamiltonian and of the applied AC voltage 
the odd and even harmonics are generated. The harmonics with largest intensity are those with 
energy close to the HOMO-LUMO gap. 
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